######### Figures in appendix ########
### Figure B16 single arima result ###

### sinlge country arimax

library(foreign)
library(TSA)
set.seed(1234567)

DE<-subset(cabinet_filter_data,countryid==1)
ES<-subset(cabinet_filter_data,countryid==2)
GB<-subset(cabinet_filter_data,countryid==3)
GR<-subset(cabinet_filter_data,countryid==4)
HU<-subset(cabinet_filter_data,countryid==5)
IE<-subset(cabinet_filter_data,countryid==6)
IT<-subset(cabinet_filter_data,countryid==7)
LV<-subset(cabinet_filter_data,countryid==8)
PL<-subset(cabinet_filter_data,countryid==9)
PT<-subset(cabinet_filter_data,countryid==10)
ROM<-subset(cabinet_filter_data,countryid==11)
AUT<-subset(cabinet_filter_data,countryid==12)
FIN<-subset(cabinet_filter_data,countryid==13)
DEN<-subset(cabinet_filter_data,countryid==14)
NET<-subset(cabinet_filter_data,countryid==15)

coef.data<-data.frame(rbind("DE","ES","GB","GR","HU","IE","IT","LV","PL","PT","ROM",
                            "AUT","FIN","DEN","NET"))
names(coef.data)<-"country"

coef.data$coef<-0
coef.data$ci.low<-0
coef.data$ci.high<-0

## North
# Denmark
mod.DEN <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
              + d.unemployment + retail+ imf+protest_freq,data=DEN) 
summary(mod.DEN)
coef.data$coef[coef.data$country=="DEN"]<- mod.DEN[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="DEN"]<- confint(mod.DEN, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="DEN"]<- confint(mod.DEN, "austerity_6m")[1,2]
box.DEN<-Box.test(mod.DEN$residuals,30,"Ljung-Box")$p.value


# Finland
mod.FIN <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
              + d.unemployment+ retail+ imf+protest_freq,data=FIN) 
summary(mod.FIN)
coef.data$coef[coef.data$country=="FIN"]<- mod.FIN[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="FIN"]<- confint(mod.FIN, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="FIN"]<- confint(mod.FIN, "austerity_6m")[1,2]
box.FIN<-Box.test(mod.FIN$residuals,30,"Ljung-Box")$p.value


# Northwest
# Austria
mod.AUT <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
              + d.unemployment+ retail+ imf+protest_freq,data=AUT) 
summary(mod.AUT)
coef.data$coef[coef.data$country=="AUT"]<- mod.AUT[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="AUT"]<- confint(mod.AUT, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="AUT"]<- confint(mod.AUT, "austerity_6m")[1,2]
box.AUT<-Box.test(mod.AUT$residuals,30,"Ljung-Box")$p.value

# Germany
mod.DE <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=DE) 
summary(mod.DE)
coef.data$coef[coef.data$country=="DE"]<- mod.DE[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="DE"]<- confint(mod.DE, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="DE"]<- confint(mod.DE, "austerity_6m")[1,2]
box.DE<-Box.test(mod.DE$residuals,30,"Ljung-Box")$p.value

# Ireland
mod.IE <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=IE) 
summary(mod.IE)
coef.data$coef[coef.data$country=="IE"]<- mod.IE[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="IE"]<- confint(mod.IE, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="IE"]<- confint(mod.IE, "austerity_6m")[1,2]
box.IE<-Box.test(mod.IE$residuals,30,"Ljung-Box")$p.value

# Netherlands
mod.NET <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
              + d.unemployment+ retail+ imf+protest_freq,data=NET) 
summary(mod.NET)
coef.data$coef[coef.data$country=="NET"]<- mod.NET[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="NET"]<- confint(mod.NET, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="NET"]<- confint(mod.NET, "austerity_6m")[1,2]
box.NET<-Box.test(mod.NET$residuals,30,"Ljung-Box")$p.value

# UK
mod.GB <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=GB) 
summary(mod.GB)
coef.data$coef[coef.data$country=="GB"]<- mod.GB[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="GB"]<- confint(mod.GB, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="GB"]<- confint(mod.GB, "austerity_6m")[1,2]
box.GB<-Box.test(mod.GB$residuals,30,"Ljung-Box")$p.value


##Southwest

# Greece
mod.GR <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=GR) 
summary(mod.GR)
coef.data$coef[coef.data$country=="GR"]<- mod.GR[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="GR"]<- confint(mod.GR, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="GR"]<- confint(mod.GR, "austerity_6m")[1,2]
box.GR<-Box.test(mod.GR$residuals,30,"Ljung-Box")$p.value


# Italy
mod.IT <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=IT) 
summary(mod.IT)
coef.data$coef[coef.data$country=="IT"]<- mod.IT[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="IT"]<- confint(mod.IT, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="IT"]<- confint(mod.IT, "austerity_6m")[1,2]
box.IT<-Box.test(mod.IT$residuals,30,"Ljung-Box")$p.value


# Portugal
mod.PT <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=PT) 
summary(mod.PT)
coef.data$coef[coef.data$country=="PT"]<- mod.PT[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="PT"]<- confint(mod.PT, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="PT"]<- confint(mod.PT, "austerity_6m")[1,2]
box.PT<-Box.test(mod.PT$residuals,30,"Ljung-Box")$p.value


# Spain
mod.ES <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=ES) 
summary(mod.ES)
coef.data$coef[coef.data$country=="ES"]<- mod.ES[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="ES"]<- confint(mod.ES, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="ES"]<- confint(mod.ES, "austerity_6m")[1,2]
box.ES<-Box.test(mod.ES$residuals,30,"Ljung-Box")$p.value


##East
# Hungary
mod.HU <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=HU) 
summary(mod.HU)
coef.data$coef[coef.data$country=="HU"]<- mod.HU[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="HU"]<- confint(mod.HU, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="HU"]<- confint(mod.HU, "austerity_6m")[1,2]
box.HU<-Box.test(mod.HU$residuals,30,"Ljung-Box")$p.value


# Latvia
mod.LV <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=LV) 
summary(mod.LV)
coef.data$coef[coef.data$country=="LV"]<- mod.LV[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="LV"]<- confint(mod.LV, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="LV"]<- confint(mod.LV, "austerity_6m")[1,2]
box.LV<-Box.test(mod.LV$residuals,30,"Ljung-Box")$p.value


# Poland
mod.PL <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
             + d.unemployment+ retail+ imf+protest_freq,data=PL) 
summary(mod.PL)
coef.data$coef[coef.data$country=="PL"]<- mod.PL[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="PL"]<- confint(mod.PL, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="PL"]<- confint(mod.PL, "austerity_6m")[1,2]
box.PL<-Box.test(mod.PL$residuals,30,"Ljung-Box")$p.value


# Romania
mod.ROM <- lm(vicab_filtered~l.vicab_filtered+austerity_6m
              + d.unemployment+ retail+ imf+protest_freq,data=ROM) 
summary(mod.ROM)
coef.data$coef[coef.data$country=="ROM"]<- mod.ROM[["coefficients"]][["austerity_6m"]]
coef.data$ci.low[coef.data$country=="ROM"]<- confint(mod.ROM, "austerity_6m")[1,1]
coef.data$ci.high[coef.data$country=="ROM"]<- confint(mod.ROM, "austerity_6m")[1,2]
box.ROM<-Box.test(mod.ROM$residuals,30,"Ljung-Box")$p.value



arima.f <- ggplot(coef.data) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_segment(aes(x=country,xend=country,y=ci.high,yend=ci.low), size=1,show.legend=F) +
  geom_point(aes(x=country,y=coef,size=2.2),show.legend=F) +
  labs(title="Country-by-country estimated immediate effect of the austerity dummy",x="",y="Effect") + 
  geom_hline(yintercept=0, linetype=2) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

arima.f

ggsave(filename="figure_B16.eps", plot=arima.f,width = 20, height = 8,device=cairo_ps)


######   Figure C1 and C2    unit root test ######

###### Check for a unit root in raw vote intention (cabinet) for each country ########## 
countrylist <- unique(country)
n <- length(countrylist)

pptest.pvalues <- rep(0,n)

for (i in 1:length(countrylist)) {
  currcty <- countrylist[i]
  
  # Check pp unit root test, omitting errors due to short series
  curpp <- try(pp.test(cabinet_filter_data$vicab_raw[country==currcty])$p.value)
  if (any(class(curpp)=="try-error")) curpp <- NA
  pptest.pvalues[i] <- curpp
}

pptest.pvalues<-data.frame(pptest.pvalues)


pptest.raw.hist<-ggplot(data=pptest.pvalues, aes(pptest.pvalues)) + 
  theme_bw()+
  geom_histogram(breaks=seq(0, 0.5, by = 0.005),binwidth = 0.01)+
  labs(title="Histogram for Phillips Perron Test p-values") +
  labs(x="p-values", y="Count")
ggsave(filename="figure_C1.eps",plot=pptest.raw.hist,width=12,height=5,device=cairo_ps)


###### Check for a unit root in filted vote intention (cabinet) for each country ########## 
countrylist <- unique(country)
n <- length(countrylist)

pptest.pvalues <- rep(0,n)

for (i in 1:length(countrylist)) {
  currcty <- countrylist[i]
  
  # Check pp unit root test, omitting errors due to short series
  curpp <- try(pp.test(cabinet_filter_data$vicab_filtered[country==currcty])$p.value)
  if (any(class(curpp)=="try-error")) curpp <- NA
  pptest.pvalues[i] <- curpp
}

pptest.pvalues<-data.frame(pptest.pvalues)

pptest.filtered.hist<-ggplot(data=pptest.pvalues, aes(pptest.pvalues)) + 
  theme_bw()+
  geom_histogram(breaks=seq(0, 0.5, by = 0.005),binwidth = 0.01)+
  labs(title="Histogram for Phillips Perron Test p-values") +
  labs(x="p-values", y="Count")
ggsave(filename="figure_C2.eps",plot=pptest.filtered.hist,width=12,height=5,device=cairo_ps)



### Figure C3 austerity_6m * dunemployment at different alphas ###
cabinet_filter_data[,"austerity_dunemployment"]<-cabinet_filter_data$austerity_6m*cabinet_filter_data$d.unemployment
attach(cabinet_filter_data)
set.seed(1234567)
int.mod.2<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_dunemployment,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.2)


## low alpha ##

int.mod.2<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_dunemployment,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.2)
alpha.low<- min(coef(int.mod.2)$l.vicab_filtered)
int.mod.2[["coefficients"]][["fixed"]][["l.vicab_filtered"]]<-alpha.low

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.25, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.5, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

Scen3 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.75, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_dunemployment<- 0
shock_scen1$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.25, na.rm = T)

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_dunemployment<- 0
shock_scen2$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.5, na.rm = T)


shock_scen3 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen3$austerity_6m [6:11]<- 1 
shock_scen3$austerity_dunemployment<- 0
shock_scen3$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.75, na.rm = T)

set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
set.seed(1234567)
Sim3 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen3, n = 24,
                   shocks = shock_scen3)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)
Sim3<-as.data.frame(Sim3)


Sim1[5,4]-Sim1[11,4]
Sim1[5,4]
Sim1[11,4]
Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=0.842",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.12, yend=-0.71), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 1st quantile of", phantom(x), Delta, "unemployment")),x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]-Sim2[11,4]
Sim2[5,4]
Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.085",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.23, yend=-0.85), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the median of", phantom(x), Delta, "unemployment")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p


Sim3[5,4]-Sim3[11,4]
Sim3[5,4]
Sim3[11,4]
Sim3.p <- ggplot(Sim3, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.398",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.37, yend=-1.02), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 3rd quantile of", phantom(x), Delta, "unemployment")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))
Sim3.p

low.alpha.dunemployment.plot<-plot_grid(NULL,Sim1.p, NULL,
                                        Sim2.p, NULL,
                                        Sim3.p,NULL, ncol=7,align = 'h',rel_widths = c(0.3, 1,0,1,0,1,0))


## high alpha ##

int.mod.2<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_dunemployment,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.2)
alpha.high<- max(coef(int.mod.2)$l.vicab_filtered)
int.mod.2[["coefficients"]][["fixed"]][["l.vicab_filtered"]]<-alpha.high

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.25, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.5, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

Scen3 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.75, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_dunemployment<- 0
shock_scen1$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.25, na.rm = T)

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_dunemployment<- 0
shock_scen2$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.5, na.rm = T)


shock_scen3 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen3$austerity_6m [6:11]<- 1 
shock_scen3$austerity_dunemployment<- 0
shock_scen3$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.75, na.rm = T)

set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
set.seed(1234567)
Sim3 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen3, n = 24,
                   shocks = shock_scen3)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)
Sim3<-as.data.frame(Sim3)


Sim1[5,4]-Sim1[11,4]
Sim1[5,4]
Sim1[11,4]
Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.279",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.18, yend=-1.09), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title="",x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]-Sim2[11,4]
Sim2[5,4]
Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.633",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.34, yend=-1.28), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title="",x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p


Sim3[5,4]-Sim3[11,4]
Sim3[5,4]
Sim3[11,4]
Sim3.p <- ggplot(Sim3, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=2.088",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.55, yend=-1.53), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title="",x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))
Sim3.p

high.alpha.dunemployment.plot<-plot_grid(NULL,Sim1.p, NULL,
                                         Sim2.p, NULL,
                                         Sim3.p,NULL, ncol=7,align = 'h',rel_widths = c(0.3, 1,0,1,0,1,0))

dunemployment.3by3<-plot_grid(low.alpha.dunemployment.plot, 
                              high.alpha.dunemployment.plot,
                              label_size = 22,
                              labels = c('Low AR (0.474)','High AR (0.697)'),
                              label_x = 0, label_y = 0, hjust = -0.05, vjust =-25.7,ncol = 1, align = 'v')

ggsave(filename="f2.eps", plot=dunemployment.3by3, width = 36, height = 13,device=cairo_ps)




### Figure C4 austerity_6m * imf at different alphas ###
cabinet_filter_data[,"austerity_imf"]<-cabinet_filter_data$austerity_6m*cabinet_filter_data$imf
attach(cabinet_filter_data)
set.seed(1234567)
int.mod.5<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_imf,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.5)


## low alpha ##

int.mod.5<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_imf,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.5)
alpha.low<- min(coef(int.mod.5)$l.vicab_filtered)
int.mod.5[["coefficients"]][["fixed"]][["l.vicab_filtered"]]<-alpha.low

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(cabinet_filter_data$protest_freq, na.rm = T),
                    austerity_imf=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(cabinet_filter_data$protest_freq, na.rm = T),
                    austerity_imf=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_imf<- 0
shock_scen1$austerity_imf[6:11]<-0

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_imf<- 0
shock_scen2$austerity_imf[6:11]<-1


set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.5, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.5, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)



Sim1[5,4]-Sim1[11,4]
Sim2[5,4]-Sim2[11,4]

Sim1[5,4]
Sim1[11,4]
Sim1[5,4]-Sim1[11,4]
Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-5,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.025",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.19, yend=-0.82), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response without external involvement")),
       x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]
Sim2[11,4]
Sim2[5,4]-Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-5,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -1.8, label = "drop=2.760",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.19, yend=-2.56), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response with external involvement")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p

low.alpha.imf.plot<-plot_grid(NULL,Sim1.p, NULL,
                              Sim2.p, NULL,
                              ncol=5,align = 'h',rel_widths = c(0.3, 1,0,1,0))


## high alpha ##

int.mod.5<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_imf,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.5)
alpha.high<- max(coef(int.mod.5)$l.vicab_filtered)
int.mod.5[["coefficients"]][["fixed"]][["l.vicab_filtered"]]<-alpha.high

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(cabinet_filter_data$protest_freq, na.rm = T),
                    austerity_imf=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(cabinet_filter_data$protest_freq, na.rm = T),
                    austerity_imf=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_imf<- 0
shock_scen1$austerity_imf[6:11]<-0

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_imf<- 0
shock_scen2$austerity_imf[6:11]<-1


set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.5, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.5, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)



Sim1[5,4]-Sim1[11,4]
Sim2[5,4]-Sim2[11,4]

Sim1[5,4]
Sim1[11,4]
Sim1[5,4]-Sim1[11,4]

Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-5,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.547",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.29, yend=-1.25), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title="",
       x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]
Sim2[11,4]
Sim2[5,4]-Sim2[11,4]

Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-5,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -1.8, label = "drop=4.229",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.29, yend=-3.93), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title="",x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p

high.alpha.imf.plot<-plot_grid(NULL,Sim1.p, NULL,
                               Sim2.p, NULL,
                               ncol=5,align = 'h',rel_widths = c(0.3, 1,0,1,0))

imf.2by2<-plot_grid(low.alpha.imf.plot, 
                    high.alpha.imf.plot,
                    label_size = 22,
                    labels = c('Low AR (0.474)','High AR (0.697)'),
                    label_x = 0, label_y = 0, hjust = -0.05, vjust =-20.7,ncol = 1, align = 'v')

ggsave(filename="figure_C4.eps", plot=imf.2by2, width = 23, height = 10,device=cairo_ps)




### Figure C5 austerity_6m * protest at different alphas ###
cabinet_filter_data[,"austerity_protest_freq"]<-cabinet_filter_data$austerity_6m*cabinet_filter_data$protest_freq
attach(cabinet_filter_data)
set.seed(1234567)
int.mod.8<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_protest_freq,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.8)


## low alpha ##

int.mod.8<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_protest_freq,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.8)

alpha.low<- min(coef(int.mod.8)$l.vicab_filtered)
int.mod.8[["coefficients"]][["fixed"]][["l.vicab_filtered"]]<-alpha.low

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.1, na.rm = T),
                    austerity_protest_freq=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.75, na.rm = T),
                    austerity_protest_freq=0)

Scen3 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.9, na.rm = T),
                    austerity_protest_freq=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_protest_freq<- 0
shock_scen1$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.1, na.rm = T)

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_protest_freq<- 0
shock_scen2$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.75, na.rm = T)


shock_scen3 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen3$austerity_6m [6:11]<- 1 
shock_scen3$austerity_protest_freq<- 0
shock_scen3$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.9, na.rm = T)

set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
set.seed(1234567)
Sim3 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen3, n = 24,
                   shocks = shock_scen3)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)
Sim3<-as.data.frame(Sim3)


Sim1[5,4]-Sim1[11,4]
Sim2[5,4]-Sim2[11,4]
Sim3[5,4]-Sim3[11,4]

Sim1[5,4]
Sim1[11,4]
Sim1[5,4]-Sim1[11,4]

Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.181",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.38, yend=-0.79), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 10% quantile of protest frequency")),x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]
Sim2[11,4]
Sim2[5,4]-Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.275",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.32, yend=-0.95), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 75% quantile of protest frequency")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p


Sim3[5,4]
Sim3[11,4]
Sim3[5,4]-Sim3[11,4]

Sim3.p <- ggplot(Sim3, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.633",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.08, yend=-1.54), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 90% quantile of protest frequency")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim3.p

low.alpha.protest.plot<-plot_grid(NULL,Sim1.p, NULL,
                                  Sim2.p, NULL,
                                  Sim3.p,NULL, ncol=7,align = 'h',rel_widths = c(0.3, 1,0,1,0,1,0))


## high alpha ##

int.mod.8<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_protest_freq,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.8)

alpha.high<- max(coef(int.mod.8)$l.vicab_filtered)
int.mod.8[["coefficients"]][["fixed"]][["l.vicab_filtered"]]<-alpha.high

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.1, na.rm = T),
                    austerity_protest_freq=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.75, na.rm = T),
                    austerity_protest_freq=0)

Scen3 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.9, na.rm = T),
                    austerity_protest_freq=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_protest_freq<- 0
shock_scen1$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.1, na.rm = T)

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_protest_freq<- 0
shock_scen2$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.75, na.rm = T)


shock_scen3 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen3$austerity_6m [6:11]<- 1 
shock_scen3$austerity_protest_freq<- 0
shock_scen3$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.9, na.rm = T)

set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
set.seed(1234567)
Sim3 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen3, n = 24,
                   shocks = shock_scen3)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)
Sim3<-as.data.frame(Sim3)


Sim1[5,4]-Sim1[11,4]
Sim2[5,4]-Sim2[11,4]
Sim3[5,4]-Sim3[11,4]

Sim1[5,4]
Sim1[11,4]
Sim1[5,4]-Sim1[11,4]

Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.698",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.55, yend=-1.13), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 10% quantile of protest frequency")),x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]
Sim2[11,4]
Sim2[5,4]-Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.853",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.46, yend=-1.38), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 75% quantile of protest frequency")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p


Sim3[5,4]
Sim3[11,4]
Sim3[5,4]-Sim3[11,4]

Sim3.p <- ggplot(Sim3, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=2.443",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.12, yend=-2.31), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 90% quantile of protest frequency")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim3.p

high.alpha.protest.plot<-plot_grid(NULL,Sim1.p, NULL,
                                   Sim2.p, NULL,
                                   Sim3.p,NULL, ncol=7,align = 'h',rel_widths = c(0.3, 1,0,1,0,1,0))

protest.3by3<-plot_grid(low.alpha.protest.plot, 
                        high.alpha.protest.plot,
                        label_size = 22,
                        labels = c('Low AR (0.474)','High AR (0.697)'),
                        label_x = 0, label_y = 0, hjust = -0.05, vjust =-25.7,ncol = 1, align = 'v')

ggsave(filename="figure_C5.eps", plot=protest.3by3, width = 36, height = 13,device=cairo_ps)